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Abstract 

n- heptane /air  premixed  turbulent  flames  in  the  high-Karlovitz  portion  of  the  thin  reaction  zone  regime  are  characterized 
and  modeled  using  Direct  Numerical  Simulations  (DNS)  with  detailed  chemistry.  To  enable  the  present  numerical  study, 
a  new  time-integration  scheme  has  been  proposed  for  the  simulation  of  stiff  reacting  flows.  Using  this  scheme,  a  series  of 
direct  numerical  simulations  of  high  Karlovitz  number,  n-CrHi6,  turbulent  premixed  flames  have  been  performed.  It  was 
found  that  the  flame  structure  of  these  turbulent  flames  can  be  well  captured  by  one-dimensional  flames  accounting  for  the 
effective  species  Lewis  numbers.  The  reaction  zone  was  found  to  remain  thin,  yet  large  fluctuations  in  the  fuel  burning  rate 
were  identified.  Extinctions  were  observed  only  in  the  presence  of  differential  diffusion,  and  these  events  were  correlated  with 
high  curvature  regions.  A  model  to  capture  the  burning  fluctuations  was  proposed  using  a  new  flamelet  approach.  For  the 
first  time,  the  evolution  of  the  turbulence  (both  turbulent  kinetic  energy  and  enstrophy)  has  been  characterized  through  the 
flame.  Under  sufficiently  high  Karlovitz  number,  the  first  Kolmogorov’s  hypothesis  has  been  confirmed.  Finally,  the  impact 
of  various  chemical/transport  model  assumptions  on  the  evolution  of  turbulence  has  been  analyzed. 


1.  Introduction 


1.1.  Regimes  of  turbulent  premixed  flames 

Turbulent  premixed  flames  are  typically  characterized  by  the  extent  to  which  turbulence  is  “expected”  to  pene¬ 
trate/disrupt  the  flame.  The  Karlovitz  number 

Ka=^,  (1) 

Tr) 

which  is  the  ratio  of  the  flame  time  scale  rp  =  If/Sl  to  the  Kolmogorov  time  scale,  provides  such  information  [1]. 
The  larger  the  Karlovitz  number,  the  more  turbulence  is  expected  to  penetrate/disrupt  the  flame. 

As  industrial  applications  of  turbulent  premixed  (and  partially  premixed)  flames  fall  in  the  thin/broken  reaction 
zone  regimes,  understanding  how  a  flame  behaves  in  these  regimes  is  critical  [2].  Experiments  are  difficult  to  conduct 
at  high  Karlovitz  numbers  and  a  limited  number  of  them  are  available  in  the  literature  [3,4].  Equivalently,  due  to 
their  expensive  computational  costs,  very  few  Direct  Numerical  Simulations  (DNS)  of  turbulent  premixed  flames  in 
the  broken  reaction  zone  or  the  thin  reaction  zone  regimes  have  been  performed  [5-10].  Figure  1  presents  all  of  these 
simulations  in  the  context  of  the  regime  diagram  (as  proposed  by  Peters  [1]).  Note  that  only  simulations  performed 
with  detailed  finite-rate  chemistry  are  presented. 

1.2.  Fuel  chemistry /transport 

The  only  fuels  that  have  been  considered  in  previous  simulations  are  hydrogen  [5,8],  methane  [6,7],  and  propane  [6]. 
While  these  fuels  are  used  in  several  ground-based  applications  [11],  most  of  the  fuels  used  for  transportation  contain 
larger  hydrocarbons  [12,13].  Since  their  chemical  pathways  are  far  more  complex,  and  a  wide  range  of  stable  species 
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Fig.  1.  DNS  of  high  Karlovitz  turbulent  premixed  flames  with  finite-rate  chemistry. 

are  present  through  the  flame  front,  it  remains  unclear  how  turbulence  influences  their  chemistry  at  high  Karlovitz 
number. 

Moreover,  heavy  hydrocarbons  have  large  Lewis  numbers  (e.g.  Le c12h26  ~  3.5).  It  has  also  been  observed  that  for 
sufficiently  high  Karlovitz  number  differential  diffusion  effects  were  negligible  for  both  H2  (Le h2  ~  0.3)  and  C3H8 
(I/ec3H8  ~  2)  [6].  However,  this  similar  behavior  between  smaller  and  larger  than  unity  Lewis  number  fuels  cannot 
be  generalized  to  intermediate  Karlovitz  numbers  as  pertinent  to  the  transition  between  the  thin/broken  reaction 
zone  regimes.  Furthermore,  while  the  series  of  lean  H2/air  premixed  flames  performed  by  Aspden  et  al.  [5]  have 
provided  information  on  how  turbulence  affects  differential  diffusion  over  a  wide  range  of  Karlovitz  numbers,  there 
is  no  such  information  available  for  heavy  hydrocarbon  fuels. 

1.3.  Approach 

To  tackle  these  questions,  a  series  of  DNS  of  a  premixed  n- C7H16  turbulent  flames  in  the  thin  reaction  zone 
regime  and  close  to  the  transition  between  the  thin/broken  reaction  zones  regimes  are  targetted  in  this  work,  n- 
C7H16  is  chosen  for  the  fuel  because  it  is  used  in  surrogates  for  gasoline  [14,13].  In  addition,  it  corresponds  to  a  step 
towards  the  simulation  of  n-decane  or  n-dodecane  flames  (larger  chemical  mechanism),  which  species  are  often  used 
in  surrogates  for  kerosene  [15,16]. 

More  specifically,  the  following  three  aspects  need  to  be  characterized  with  these  target  flames: 

(i)  the  effect  of  turbulence  on  the  flame  structure  and  the  impact  of  turbulent  mixing  on  differential  diffusion, 

(ii)  the  effect  of  turbulence  on  the  reaction  zone  with  and  without  differential  diffusion, 

(iii)  the  effect  of  the  flame  on  the  evolution  of  the  turbulent  flow  field. 

2.  Governing  equations 

2.1.  Fluid  mechanics 

The  reacting  mixture  is  assumed  to  contain  a  total  number  of  N  species  and  their  chemistry  is  assumed  to  be  given 
by  a  chemical  kinetics  mechanism  involving  K  reactions,  with  forward  and  backward  reactions  counted  separately. 
The  chemically  reacting  flows  of  interest  in  the  current  study  are  of  relatively  low  Mach  number  (Ma),  typically  below 
0.3  [17,18].  Under  this  condition,  the  acoustic  waves  can  be  ignored  and  the  pressure  field  can  be  decomposed  into 
a  spatially-invariant,  but  (potentially)  time-dependent  component,  Po  (£),  and  a  fluctuating  hydrodynamic  pressure, 
p  (x,  t)  [17-19,1].  Since  the  focus  is  on  turbulence-chemistry  interaction,  Soret  and  Dufour  effects,  body  forces,  and 
radiative  heat  transfer  are  ignored  [20,21,18,22]. 

Under  these  assumptions,  the  evolution  of  the  system  is  governed  by  the  following  conservation  equations  of  mass, 
momentum,  energy,  and  species  density  [23,19,1]: 

%  +  V  •  (pu)  =  0  (2) 

dpu  „  .  .  „  „ 

——  +  V  •  (pu  0  u)  =  -  Vp  +  V  •  r 
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(3) 


(4) 


dpT 

dt 


+  V  •  (/ouT) 


V  •  (pcpaVT)  +  ^  cP)ip  (-^- VF;  +  YSVc,!  )  •  VT 


'  +  Co T 


®  +  V  •  (puF*)  =  V  •  (p^VFij  +  V  •  (pFiVci)  +  Wi.  (5) 

In  the  above  equations,  p  is  the  density,  u  is  the  velocity  vector,  T  denotes  the  temperature  of  the  mixture,  and  Yj 
is  the  mass  fraction  of  species  i.  In  the  momentum  equation  (Eq.  3),  r  is  the  deviatoric  stress  tensor,  defined  as 


t  =  M[Vu+(Vu)t]  -^(V-u)I, 


(6) 


where  I  is  the  identity  matrix  and  p  is  the  fluid  viscosity.  In  the  energy  conservation  equation  (Eq.  4),  lot  includes 
heat  source  terms  due  to  chemical  reactions,  a  is  the  thermal  diffusivity,  and  cp  is  the  specific  heat  at  constant 
pressure  of  the  mixture,  given  by 

N 

Cp  =  ^  ^  YiCpti,  (7) 

where  cp ^  is  the  specific  heat  at  constant  pressure  of  species  i.  In  the  species  conservation  equations  (Eq.  5),  LOj  is 
the  chemical  source  term  of  species  i,  and  Lej  is  the  Lewis  number  of  species  i,  defined  as 


Lei  =  t 

with  Di  the  mass  diffusivity  for  species  i.  The  correction  velocity  Vc  ;  in  Eq.  5  accounts  for  gradients  in  the  mixture 
molecular  weight  as  well  as  ensures  zero  net  diffusion  flux.  It  has  the  following  expression  [19,1]: 


VC)i  = 


a  VW 

Yei~W~ 


N 


-HE 


where 


VYh 


•  i  -^Cj 
3=1 


N 


—  a- 


VW 

w~ 


N 


Y, 


<3  =  1 


-1 


w=  lyL 

l  ^  Wj 


K3=l 


is  the  local  mean  molecular  weight  of  the  mixture,  and  Wj  is  the  molecular  weight  of  species  j. 
The  above  set  of  equations  is  complemented  by  the  equation  of  thermodynamic  state 


P  = 


PqW 

RT  ’ 


where  Pq  is  the  thermodynamic  pressure  and  R  is  the  universal  gas  constant. 


(9) 


(10) 


(11) 


2.2.  Chemical  model 


The  overall  rate  of  change  of  species  i,  cH,  in  Eq.  5  can  be  split  into  a  production  term,  c of ,  and  a  consumption 
term,  c jr,  as 

LOi  =  Cof  —  .  (12) 

The  production  rate  of  species  i,  Cof is  given  by  the  sum  of  the  contributions  from  all  elementary  chemical  reactions 
leading  to  the  formation  of  this  species: 


Cjf  =  Wj 


K 

E 

3  =  1 

V  j  i  0 


N 


k>  n 


(13) 


where  Vjs  is  the  stoichiometric  coefficient  of  species  s  in  reaction  j.  In  the  above  expression,  the  rate  constant  of 
reaction  j,  kj ,  is  given  by  the  Arrhenius  form,  kj  (T)  =  AjTbi  exp_Ta^/T,  where  Taj  is  the  activation  temperature 
of  this  reaction.  The  consumption  rate  is  writen  in  a  similar  manner. 

The  local  heat  release  rate  is  given  by 

N 

*r*<*5>*«’  (14) 

3= 1 


where 
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(15) 


hi  —  Y$  + 


f 

J  Tq 


CpidT 


is  the  specific  enthalpy  of  species  i,  and  h 3  denotes  its  value  under  standard  and  reference  conditions. 


3.  Numerical  algorithm 

3.1.  Overview  of  the  numerical  solver 

The  simulations  in  this  work  are  performed  using  the  structured,  multi-physics  and  multi-scale  finite- difference 
code  NGA  [17].  The  NGA  code  allows  for  accurate,  robust,  and  flexible  simulations  of  both  laminar  and  turbulent 
reactive  flows  in  complex  geometries  and  has  been  applied  in  a  wide  range  of  test  problems,  including  laminar 
and  turbulent  flows  [23-25]  and  constant  and  variable  density  flows  [17,26,27],  as  well  as  Large-Eddy  Simulations 
(LES)  [24,28]  and  Direct  Numerical  Simulations  (DNS)  [29,27,30].  This  numerical  solver  has  been  shown  to  conserve 
discretely  mass,  momentum,  and  kinetic  energy,  with  arbitrarily  high  order  spatial  discretization  [17].  NGA  uses 
both  spatially  and  temporally  staggered  variables  [17].  All  scalar  quantities  (p,p,  T,  D{ ,  a,  cp,  cPj,  /a)  are  stored  at 
the  volume  centers,  and  the  velocity  components  are  stored  at  their  respective  volume  faces.  The  convective  term  in 
the  species  transport  equations  is  discretized  using  the  bounded  quadratic  upwind-biased  interpolative  convective 
scheme  (BQUICK)  [31],  and  the  diffusive  term  is  discretized  using  a  second-order  centered  scheme.  The  variables 
are  advanced  in  time  using  the  second-order  semi-implicit  Crank-Nicolson  scheme  of  Pierce  and  Moin  [32]. 

3.2.  Iterative  procedure 


An  iterative  procedure  is  applied  to  fully  cover  the  non-linearities  in  the  Navier-Stokes  equations.  This  iterative 
procedure  has  been  found  to  be  of  critical  importance  for  stability  and  accuracy  considerations  [17,18,32].  The 
density,  pressure,  and  scalar  fields  are  advanced  from  time  level  £n+V2  to  £n+3/2,  and  the  velocity  fields  are  advanced 
from  time  level  tn  to  £n+1.  The  subscript  k  refers  to  the  sub-iterations. 

Traditionally,  the  scalar  fields  are  advanced  in  time  using  the  semi-implicit  Crank-Nicolson  method  [17,32]  for  the 
convective  and  diffusive  terms,  and  explicit  integration  for  the  chemical  source  terms: 


(16) 


7i+3/2-y-7i+3/2 

Pk  1  k- fs| 


=  /)n+l/2Yn+l/2  +  ^ 


[(C 


n- jMt 
k 


■«2 


+ 


At 

~2 


fdC  <9D\n+1 
(3Y  +  3Y,)fc 


f  -*^-77+3/2  -*^-77+3/2 

\  1  fc+1  —  1  k 


(17) 


To  simplify  the  discrete  notations  for  spatial  differential  operators,  the  operators  corresponding  to  the  convective 
and  diffusive  terms  in  the  scalar  equations  (Eq.  5)  are  written  as  C  and  D,  respectively.  an(J  §Y  are  Jacobian 
matrices  corresponding  to  the  convective  and  diffusive  terms,  respectively.  For  simpler  implementation,  the  set  of 
equations  (Eq.  17)  is  solved  in  practice  in  its  residual  form: 


y«+3/2  =  Yn+3/2  _  A/  J  1 


*  ©fc, 


(18) 


where  the  matrix  J  is  defined  as 
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and  the  vector 
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(19) 

(20) 


is  the  error  (residual)  made  on  the  species  transport  equation  at  the  previous  sub-iteration.  When  the  sub-iterations 
are  fully- converged,  the  residual,  ©&,  is  zero. 


3.3.  Preconditioned  iterative  method 


The  time- marching  for  species  transport  equations  described  in  Eq.  18  resembles  the  standard  preconditioned 
Richardson-type  iterative  method  [33],  where  the  matrix  J  acts  as  a  preconditioner.  More  precisely,  the  choice  of  the 
preconditioner,  J,  can  be  arbitrary  and  does  not  modify  the  discrete  form  of  the  equations  to  solve  (i.e.  =  0).  It 

only  changes  the  convergence  characteristics  of  the  iterative  method.  For  instance,  setting 
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J  =  pnk+3/2I  (21) 

is  equivalent  to  the  fully-explicit  integration  of  the  convective,  diffusive,  and  chemical  source  terms  in  the  species 
transport  equations,  while  setting 


,  _  n+3/2T  At  /  dc  dT>  on  \ n+1 
3~Pk  I~Y\dY  +  dY  +  dY)k 


(22) 


corresponds  to  the  full-implicit  integration  of  the  convective,  diffusive,  and  chemical  source  terms  in  the  species 
transport  equations. 

Clearly,  there  is  a  trade-off  in  the  choice  of  the  preconditioner.  Since  it  is  applied  at  each  step  of  the  iterative 
method,  it  is  preferable  to  have  a  preconditioning  matrix,  J,  with  low  computing  and  inversion  cost.  The  cheapest 
preconditioner  would  therefore  be  the  one  described  by  Eq.  21  (fully-explicit  integration),  which  leads  to  poor 
convergence  performance  requiring  extremely  small  time  step  sizes.  On  the  other  hand,  the  optimal  preconditioner 
would  be  the  one  leading  to  the  fully-implicit  integration  of  the  various  terms  (Eq.  22).  Unfortunately,  since  the 
chemical  source  terms  of  most  species  are  generally  dependent  on  a  large  number  of  other  species,  the  chemical 
Jacobian  matrix,  is  usually  not  sparse  [34].  Therefore,  its  construction  and  inversion  may  become  prohibitively 
expensive  especially  when  a  large  number  of  species  is  considered. 

The  proposed  preconditioner  is  written  as 


j  =  Pnk+3/2i 


At  fdC  dB  _  \n+1 

TUy  +  5Y“AT  ’ 


(23) 


where  A  is  a  diagonal  matrix  defined  as 


(24) 


The  matrix  A  may  be  regarded  as  a  very  good  approximation  of  the  diagonal  of  the  chemical  Jacobian  (see  Fig.  2). 
The  proposed  preconditioner  aims  to  suppress  the  small  timescales  due  to  the  fast  consumption  of  the  different 
species  in  the  system  with  stiff  chemistry. 
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Fig.  2.  Comparison  of  the  chemical  timescale  (r)  of  the  full  chemical  Jacobian  to  the  species  lifetime  of  the  preconditioned  chemical 
Jacobian  at  the  peak  rate  of  heat  release,  considering  various  fuels  and  chemical  mechanisms  and  temporal  accuracy  of  the  method. 


3.4.  Scheme  performance 

The  performance  of  the  proposed  method  was  numerically  tested  on  two  flow  configurations:  a  one-dimensional 
unstretched  premixed  flame  and  a  three-dimensional  turbulent  premixed  flame,  both  with  an  unburnt  mixture  of  air 
and  n- heptane.  The  method  was  shown  to  be  second-order  accurate  in  time  (Fig.  2).  The  computational  cost  of  a 
single  iteration  with  the  proposed  method  is  similar  to  that  of  an  explicit  time-integration  scheme.  Therefore,  the 
simulation  speed-up  achieved  with  the  proposed  method  corresponds  to  the  increase  in  the  largest  stable  time  step 
size.  For  the  three-dimensional  turbulent  premixed  flame,  the  simulation  could  be  performed  with  a  convective  CFL 
of  0.8  (optimal,  with  or  without  chemistry). 

The  theoretical  analysis  for  stability  and  convergence  rate  is  general  and  is  not  limited  by  the  type  of  fuel,  chemical 
mechanism,  or  flow  configuration.  Therefore,  it  was  repeated,  in  the  context  of  one-dimensional  premixed  flames, 
with  several  fuels,  unburnt  conditions,  and  chemical  mechanisms.  It  was  also  performed  with  non-premixed  flamelets 
using  different  scalar  dissipation  rates.  The  method  provides  good  convergence  rates  of  the  sub- iterations  close  to 
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Fig.  3.  Computational  domain  demonstrating  the  approximate  location  of  the  flame  and  region  of  forcing.  Diagram  taken  from  Ref.  [38]. 


A 

B 

c 

C' 

D 

E 

Tu  [K] 

298 

298 

800 

500 

800 

298 

Put  Pb 

7.8 

7.8 

3.3 

4.9 

3.3 

7.8 

Lx[  mm] 

25.6 

25.6 

16.8 

18.7 

16.8 

25.6 

L  [mm] 

2.33 

2.33 

1.53 

1.70 

1.53 

2.33 

Grid 

11  x  1283 

11  x  1283 

11  x  1283 

11  x  1463 

11  x  2203 

11  x  2403 

Ret 

83 

190 

170 

290 

380 

390 

K  CLU 

70 

220 

200 

650 

750 

640 

u'/SL 

9 

18 

19 

38 

45 

37 

Io/If 

1.1 

1.1 

1.2 

1 

1.2 

1.1 

PuIpM 

16 

9 

7 

4.6 

3.5 

5.1 

Sl  [m/s] 

0.36 

0.36 

2.3 

0.86 

2.3 

0.36 

l  p  [mm] 

0.39 

0.39 

0.25 

0.32 

0.25 

0.39 

Table  1 

Physical  and  numerical  parameters  of  the  DNS  which  employ  finite-rate  chemistry  and  non-unity  Lewis  number  transport. 

the  stability  limit  for  all  the  chemical  mechanisms  considered.  The  proposed  preconditioning  method  showed  great 
potential  for  the  efficient  time-integration  of  turbulent  flames.  Although  not  a  primary  target,  the  method  was  also 
shown  to  work  for  a  homogeneous  ignition  case. 

Finally,  the  proposed  method  is  more  suited  than  other  methods  for  reacting  flows  in  which  the  convective 
timescales  are  of  the  order  of  10-6  s  or  less.  These  correspond  to  moderately  to  highly  turbulent  (non-premixed  or 
premixed)  flames  (high  Karlovitz  for  premixed  flames). 

4.  Direct  Numerical  Simulations 


qCc0  9GcAO 


o  o  o  o  o 
O  O  Q  o 


o 


outflow 


4.1.  Physical  Configuration 

The  present  study  considers  statistically-stationary,  statistically-planar  premixed  turbulent  n-heptane/air  flames 
at  a  slightly  lean  equivalence  ratio  =  0.9)  and  atmospheric  pressure.  The  three-dimensional  domain  has  an 
inflow  and  outflow  at  the  left  and  right  x  boundaries,  respectively,  and  periodic  boundary  conditions  in  the  y  and  z 
directions  (Fig.  3).  The  height  and  width  of  the  channel  are  equal  and  denoted  as  L,  while  the  length,  Lx ,  is  equal 
to  11 L.  A  separate  DNS  is  performed  of  relatively  weak,  homogeneous,  isotropic,  triply  periodic  box  turbulence  and 
is  used  to  generate  the  inflow  condition.  The  mean  inflow  velocity  is  constant  for  each  case  and  set  to  a  value  which 
approximates  the  turbulent  flame  speed,  allowing  for  an  arbitrary  long  run-time.  This  configuration  lacks  any  mean 
shear  so  that  the  effects  of  the  flame  on  the  turbulence  may  be  specifically  studied. 

The  turbulence  and  temperature  in  the  reactants  are  varied  between  simulations  to  investigate  the  effects  of  both 
the  unburnt  Karlovitz  number  and  flame  density  ratio.  All  necessary  information  about  the  different  simulations 
is  provided  in  Table  1  and  2,  where  Ret  =  u'l0/v ,  u! ,  and  r]u  are  the  turbulent  Reynolds  number,  rms  velocity 
fluctuation,  and  Kolmogorov  length  scale,  respectively,  all  calculated  in  the  unburnt  gas.  Cases  Bi,  Bxab,i,  Bos,r 
and  B^ab  x  are  performed  to  test  the  effects  of  the  transport  models,  chemical  models,  and  Reynolds  number,  and 
use  the  same  method. 

Between  the  cases  studied,  the  unburnt  Karlovitz  number  varies  by  an  order  of  magnitude  (Kau= 70  -  760)  and 
the  unburnt  temperature  spans  all  practically  relevant  conditions  (Tu= 298  -  8007T).  Cases  A,  B,  and  E  have  the 
same  density  ratio  and  are  used  to  test  the  effects  of  Kau  independently,  while  the  pairs  B,  C  and  E,  C'  have  the 
same  Kau  and  are  used  to  test  the  density  ratio  independently.  Figure  4  shows  Kau  and  the  density  ratio  for  each 
case  as  well  as  their  location  on  the  Peters’  regime  diagram.  These  conditions  span  the  transition  from  the  thin  to 
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Bi 

Bxab,l 

Bos,i 

r)4 

°Tab,l 

Tu  [K] 

298 

298 

298 

298 

Pu/  Pb 

7.8 

7.8 

7.3 

7.8 

Lx[  mm] 

25.6 

25.6 

25.6 

60.6 

L  [mm] 

2.33 

2.33 

2.33 

9.32 

Grid 

11  x  1283 

11  x  1283 

11  x  1283  2574  x  5122 

Ret 

190 

190 

190 

1150 

Kau 

280 

280 

250 

280 

u'/SL 

21 

21 

22 

33 

IJIf 

1 

1 

1.2 

4 

V u  [^m] 

9 

9 

9 

9 

Sl  [m/s] 

0.29 

0.29 

0.27 

0.29 

Mm] 

0.43 

0.43 

0.36 

0.43 

Table  2 

Physical  and  numerical  parameters  of  additional  DNS  which  vary  the  transport  and  chemical  models.  Subscripts  1,  Tab,  and  OS  corre¬ 
spond  to  simulations  using  unity  Lewis  numbers,  tabulated  chemistry,  and  one-step  chemistry  respectively.  Superscript  of  4  corresponds 
to  Iq/If  —  4. 


3 


Pu/  Pb 

(a)  Karlovitz  numbers  and  density  ratios 


l /If 

(b)  Peters  regime  diagram 


Fig.  4.  Conditions  in  the  unburnt  flow  of  performed  simulation.  Symbols  on  the  left  plot  correspond  to  the  same  simulations  on  the  right. 


broken/distributed  reaction  zone  regimes. 

Each  simulation  is  run  for  over  25  eddy  turnover  times  (r0  =  fc/e,  where  k  is  the  TKE)  after  an  initial  transient 
period,  in  order  to  provide  sufficient  statistical  samples.  Further  specifications  of  the  simulation  conditions  are  listed 
in  Table  1  and  2. 


4.2.  Turbulence  forcing 

As  all  other  similarly  high  Ka  numerical  simulations  from  the  literature  [5,6,39],  the  present  configuration  misses 
the  generation  of  turbulence  due  to  large  scale  flow  straining  (larger  than  the  domain  size).  As  a  result,  the  turbulence 
is  expected  to  decay,  and  the  use  of  velocity  field  forcing  is  necessary. 

In  previous  work,  spectral  forcing  techniques  were  often  used  to  offset  the  decay  of  TKE  and  maintain  the 
turbulence  characteristics  [5,39].  In  the  present  work,  the  linear  velocity  forcing  method  [40,41,30]  was  preferred  for 
its  more  physical  nature  and  good  stability  properties  [26].  The  linearly  forced  turbulent  field  under  comparable 
Reynolds  number  was  analyzed  in  Ref.  [30]  and  it  was  shown  that  the  second-  and  third-order  structure  functions 
and  the  energy  spectrum  are  self-consistent  and  in  agreement  with  experimentally  obtained  data  [42]  of  decaying 
grid  turbulence.  The  linear  forcing  method  mimics  the  missing  large  scale  straining  by  appending  a  source  term  to 
the  momentum  equation  [40,41]. 

Figure  5  presents  the  energy  spectra  in  the  unburnt  gases  for  both  the  lower-  and  the  higher-ITa  flames.  First, 
it  is  important  to  note  that,  as  expected,  both  profiles  collapse  on  each  other.  Second,  the  presence  of  an  inertial 
sub-range  is  very  limited.  This  is  to  be  expected  given  the  turbulent  Reynolds  numbers.  Finally,  for  both  cases,  the 
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turbulent  kinetic  energy  is  contained  over  two  decades  of  length  scales.  This  is  also  to  be  expected  given  the  turbulent 
Reynolds  numbers  of  approximately  100.  This  range  of  length  scale  is  limited  by  the  inherent  computational  cost. 


KT| 


Fig.  5.  Normalized  two-dimensional  three  components  spectra  taken  in  a  y-z  plane  in  the  unburnt  gases  (averaged  over  time).  K  =  2kir  /  L , 
for  k  =  1,  2, ...,  Ny/ 2,  with  Ny  the  number  of  points  in  the  y-  or  ^-direction,  is  the  wavenumber. 


4.3.  Forcing  analysis 


As  discussed  previously,  the  implemented  method  of  turbulence  forcing  intends  to  maintain  a  constant  TKE 
throughout  domain.  Figure  6a  displays  the  planar  averaged  TKE  for  case  B,  which,  as  expected,  is  nearly  constant 
in  the  region  of  forcing.  Next  we  consider  the  planar  averaged  dissipation  rate,  defined  as  e  =  r'  :  S" /p  [43].  Figure  6b 
shows  that  this  quantity  is  also  constant  in  the  region  of  forcing.  This  may  be  explained  through  the  TKE  transport 
equation  for  the  case  of  statistical  stationarity  and  homogeneity  in  the  y  and  z  directions, 


u—  =  —  e  +  2  Ak0  —  u"u" 
ox 


du 

dx 


1  d 
+  p  dx 


-u"p(u"  •  u")  -  P'u' 


i  _  u"  BP  1 _ 

-t  :  S"  -  —  4-  +  =P'V  •  u" 

P  p  OX  p 


(25) 

(26) 


where  ex  is  the  unit  vector  in  the  x  direction.  The  first  four  terms  in  this  equation,  namely  the  left  hand  side  (LHS), 
dissipation,  forcing,  and  dilatation,  respectively,  are  plotted  in  Fig.  6c  along  with  with  the  residual,  which  represents 
the  cumulative  magnitude  of  all  the  remaining  terms.  Dissipation  and  forcing  are  the  dominant  terms,  even  through 
the  flame. 

Previous  studies  in  3D  periodic  homogeneous,  isotropic  turbulence  found  that  linear  forcing  results  in  an  integral 
length  scale,  l  =  (2&/3)3/2/e,  which  is  proportional  to  the  domain  size  [44],  namely  about  0.19L.  The  integral  length 
scale  for  the  present  configuration  is  fairly  constant  through  the  domain  and  acquires  a  value  of  approximately  0.16L 
(Fig.  6d). 


4.4.  Case  B  -  Results  overview 

At  high  Karlovitz  numbers,  the  spatial  structure  of  a  premixed  flame  is  expected  to  depart  from  that  of  a  laminar 
flame,  as  turbulent  structures  are  sufficiently  small  enough  to  penetrate  the  preheat  zone,  and  potentially  the 
reaction  zone.  Figure  7  presents  contours  of  vorticity,  n-C7Hi6  and  CH2O  mass  fractions,  and  temperature  through 
the  flame  brush  on  a  two-dimensional  horizontal  slice  of  the  non-unity  Lewis  number  turbulent  flame.  This  figure  is 
representative  of  both  the  unity  and  the  non-unity  Lewis  number  turbulent  flame  brushes. 

As  expected  the  preheat  zone  appears  largely  thickened  [1]  (approximately  10  times  larger  than  the  laminar  flame). 
It  is  interesting  to  note  that  smaller  turbulent  structures  are  observed  upstream  of  the  flame  (i.e.  in  the  preheat 
zone)  compared  to  close  to  the  reaction  zone.  This  is  partially-explained  by  the  fact  that  the  kinematic  viscosity 

increases  (by  up  to  a  factor  of  30)  and  the  Kolmogorov  length  scale,  77  =  (^/e)1^4,  increases  through  the  flame  by 
about  a  factor  of  13. 
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(a)  Turbulent  kinetic  energy.  (b)  Dissipation  rate. 


0  10  20  30  40  50  60 

x  Jlp 

(c)  TKE  budget. 


(d)  Integral  length  scale. 


Fig.  6.  Planar  and  temporal  average  of  turbulent  kinetic  energy,  dissipation  rate,  and  terms  in  the  TKE  transport  equation,  and  integral 
length  scale  for  case  B.  Green  dashed  lines  correspond  to  averages  when  either  the  first  or  second  half  of  the  data  is  used;  they  are 
indicative  of  the  statistical  uncertainty  in  the  computed  averages. 


5.  Effective  Lewis  numbers 

5.1.  Turbulent  flame  structure  in  the  absence  of  differential  diffusion 

To  properly  assess  the  influence  of  turbulence  on  the  flame  structure,  one  can  analyze  the  correlation  between 
species  and  temperature  (or  any  other  progress  variable).  As  such,  the  flame  structure  can  be  adequately  compared 
to  that  of  a  one-dimensional  laminar  flame,  which  is  well  represented  in  temperature  space.  Any  departure  from  a 
one-dimensional  laminar  flame  due  to  turbulence  should  be  captured  by  these  species  mass  fraction  profiles. 

In  this  sense,  several  species  mass  fractions  are  plotted  against  temperature  and  are  compared  to  their  one¬ 
dimensional  laminar  flame  equivalent.  Figure  8  shows  joint  probability  densities  of  n-CyHie,  C2H4,  and  CO2  mass 
fraction,  vs.  temperature.  These  species  correspond  to  a  reactant,  an  intermediate  species,  and  a  product,  respectively. 
The  conditional  mean  of  these  species  mass  fraction  (conditional  on  temperature)  is  also  shown.  This  figure  is 
representative  of  the  overall  flame  structure  as  the  mass  fractions  of  other  species  show  similar  behaviors.  These 
results  suggest  that  the  influence  of  turbulence  on  the  flame  structure  in  the  absence  of  differential  diffusion  is 
very  limited  as  the  spread  of  the  joint  PDF  is  limited  (this  has  also  been  observed  by  Aspden  et  al.  for  a  high- Ka 
CHU/air  flame  [6]).  More  interestingly,  the  conditional  mean  profiles  of  these  species  follow  very  closely  the  profiles 
of  a  one-dimensional,  unstretched  laminar  flame  at  the  same  condition.  This  result  is  surprising,  as  the  turbulent 
flame  is  clearly  not  in  the  flamelet  regime  and  does  not  look  like  a  flamelet. 

In  summary,  while  the  turbulent  flame  is  thickened  by  turbulence  and  is  clearly  not  a  thin  flame,  its  structure  is 
similar  to  that  of  a  flamelet.  This  may  suggest  that  the  use  of  a  progress  variable  with  tabulated  chemistry  [45-47] 
would  be  justified  and  sufficient  even  at  such  high  Karlovitz  number. 
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Fig.  7.  Contours  of  vorticity,  n-CjYiiQ  and  CH2O  mass  fractions,  and  temperature  through  the  flame  brush  on  a  two-dimensional 
horizontal  slice.  The  laminar  flame  thickness  Ip  is  added  for  comparison. 
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Fig.  8.  Joint  PDF  and  conditional  mean  (solid  line)  of  the  71-C7H16  (top  left),  C2H4  (top  right),  and  CO2  (bottom)  mass  fraction  vs. 
temperature  from  the  unity  Lewis  number  DNS.  The  unity  Lewis  number  flamelet  solution  is  also  shown  (dashed  line). 


5.2.  Turbulent  flame  structure  with  differential  diffusion 


Similarly  to  Fig.  8,  Fig.  9  presents  the  structure  of  C2H4  through  the  non-unity  Lewis  number  flame.  The  full- 
transport  flamelet  solution  is  also  added  for  comparison.  While  turbulence  has  almost  no  impact  on  the  structure 
of  the  unity  Lewis  number  flame,  it  has  a  clear  effect  on  that  of  a  non-unity  Lewis  numbers  flame.  In  this  case,  the 
turbulent  flame  structure  lies  between  that  of  a  full  transport  and  a  unity  Lewis  number  flamelet.  Once  again,  a 
first  order  effect  of  turbulence  on  scalar  transport  is  an  increase  in  the  effective  diffusivity  of  each  scalar  (including 
species  mass  fractions  and  temperature)  through  increased  mixing.  As  a  result,  the  effective  species  Lewis  number 
take  the  following  from: 


Leie  ff 


ot  -j-  Dp 
Di  +  Dt 


(27) 


A  similar  expression  for  these  effective  Lewis  numbers  was  first  suggested  by  Peters  [1].  Equation  27  suggests  that  if 
the  turbulence  were  sufficiently  intense,  the  non-unity  Lewis  number  case  would  behave  the  same  as  the  unity  Lewis 
number  case,  i.e.  turbulence  would  suppress  differential  diffusion  effects. 


5.2.1.  Effective  Lewis  numbers 

Any  model  should  at  least  capture  the  mean  of  the  turbulent  flame  (in  a  statistical  sense).  Reynolds-averaged 
transport  equations  provide  information  about  the  ensemble/statistical  average  of  the  transported  quantity: 
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Fig.  9.  Joint  PDF  and  conditional  mean  (solid  line)  of  the  C2H4  mass  fraction  vs.  temperature  from  the  non-unity  Lewis  number  DNS. 
The  non-unity  and  unity  Lewis  number  flamelet  solutions  are  also  shown  (dashed  line). 


d_ 

dt 


+  V- 


=  V 


p(Di  +  DT)VYi 


(28) 


p  (<a  +  Dt)  VT 


(29) 


where  7  denotes  a  Reynolds  average  and  r  a  Favre  average.  DT  is  the  eddy  diffusivity  and  is  assumed  to  be 
identical  for  all  species  and  temperature.  This  result  is  a  direct  consequence  of  assuming  turbulence  mixes  all  scalars 
the  same  way.  A  similar  assumption  is  made  in  transported  PDF  methods  [48,49].  From  Eq.  28  and  29,  an  effective 
Lewis  number  can  be  obtained  for  each  species: 


Di 


heat 


Dt 


D„ 


1+Dt 


Di  +  DT  -T  +  —' 

1  1  La  1  a 


(30) 


Assuming  that  the  eddy  diffusivity  scales  as  the  product  of  the  integral  length  scale  and  the  velocity  fluctuation, 
Dt  oc  lur ,  the  model  for  the  effective  Lewis  numbers  can  be  expressed  as 


r  RANS 
^ez,eff 


1 


,RANS 


ReT 


^+aRANSReT’ 


(31) 


where  a  proportionality  constant  aRANS  is  used.  With  this  approach,  the  effective  Lewis  number  does  not  depend 
on  any  flame  characteristic.  This  result  is  rather  surprising  as  such  a  dependency  should  be  expected. 

Assuming  that  the  eddy  diffusivity  is  controlled  by  eddies  of  the  size  of  the  flame,  Dt  oc  Ipup,  a  model  for  the 
effective  Lewis  numbers  follows 

1  +  alFKa2/3 


= 


jD+clDKcl2/ 3’ 


(32) 


where  a}F  is  a  proportionality  constant. 

A  third  approach  is  to  consider  the  response  time  of  the  flame  to  unsteady  perturbations.  This  response  time  is 
expected  to  scale  like  the  flame  time  tp,  i.e.  Dt  oc  l(tp)u(tF).  The  model  for  the  effective  Lewis  numbers  becomes 


tF  l  +  a^Ka* 

*’6ff  Ze~  +  atF Ka2  ’ 


(33) 


where  a  proportionality  constant  atp  is  used. 

Finally,  it  was  found  that  an  additional  empirical  model  in  which  the  ReT  is  simply  replaced  by  Ka  performs 
better  for  a  series  of  hydrogen/air  turbulent  premixed  flames.  This  model  is  given  by 


Le 


Ka  _ 
i,eff  — 


1  +  aKaKa 
+  aKaKa 


(34) 


5.3.  Comparison  with  DNS  data 

The  iTa-based  model  (Eq.  34)  is  used  to  compute  a  set  of  effective  species  Lewis  numbers.  The  solution  of  an 
unstretched  laminar  flame  with  this  set  of  effective  species  Lewis  numbers  is  obtained  with  FlameMaster.  The  species 
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Fig.  10.  Conditional  mean  (solid  line)  of  the  71-C7H16  (top  left),  C2H4  (top  right),  and  CO2  (bottom)  mass  fraction  vs.  temperature. 
The  non-unity,  the  unity,  and  the  effective  Lewis  number  flamelet  solutions  are  also  shown  (dashed  line). 


mass  fraction  vs.  temperature  profiles  are  compared,  in  Fig.  10,  to  the  corresponding  conditional  means  from  the 
DNS  for  the  fuel,  C2H4,  and  CO2.  Although  differences  are  noticeable,  especially  for  the  product,  the  effective  Lewis 
number  model  captures  fairly  well  the  effect  of  turbulent  mixing  on  the  mean  flame  structure. 


6.  Reaction  zone 

6.1.  Overview 

Contours  of  the  source  term  of  n- C7H16  and  H2O  for  the  non-unity  Lewis  number  B  flame  are  shown  in  Fig.  11  on 
the  same  two-dimensional  slice  as  in  Fig.  7.  For  a  Karlovitz  numbers  as  large  as  220,  scaling  arguments  suggest  that 
turbulent  structures  should  be  sufficiently  small  enough  to  penetrate  the  reaction  zone  [1].  However,  two  observations 
can  be  made  from  Fig.  11:  1)  the  reaction  zone  appears  to  be  thin  and  2)  signs  of  local  extinctions  (broken  reaction 
zone)  can  be  observed  in  the  fuel  consumption  zone,  but  cannot  be  observed  in  the  H2O  production  zone. 

The  first  observation  could  seem  inconsistent  with  the  fact  that  the  Kolmogorov  length  scale  in  the  unburnt  gas 
is  10  times  smaller  than  the  laminar  reaction  zone  thickness.  However,  as  mentioned  above,  the  Kolmogorov  length 
scale  increases  through  the  flame  due  to  the  increasing  kinematic  viscosity,  and  becomes  as  large  as  the  laminar 
reaction  zone  thickness  as  it  reaches  the  burnt  side. 


wh2o 

(kg  m'3  s'1) 


250 


125 


Fig.  11.  Contours  of  the  source  terms  of  71-C7H16  (left)  and  H2O  (right)  on  the  same  two-dimensional  horizontal  slice  as  in  Fig.  7.  The 
isoterm  T  =  1500  K  (white)  is  also  shown.  The  laminar  reaction  zone  thicknesses  (full  width  at  half-height)  of  71-C7H16  (5c7Hi6)  an<l 
H2O  (<5h2o)  are  also  shown  for  comparison. 


6.2.  Fuel  consumption  at  Tpeak 

To  investigate  the  impact  of  turbulence  on  the  reaction  zone,  Fig.  12(a)  shows  the  joint  probability  density  of  the 
normalized  fuel  consumption  rate  as  a  function  of  temperature  for  the  unity  Lewis  number  flame.  The  mean  of  the 
source  term  conditional  to  temperature  is  also  shown.  This  profile  is  compared  to  its  unstretched  one-dimensional 
flame  equivalent.  While  the  mean  profile  is  very  close  to  that  of  the  one-dimensional  flame,  fluctuations  around  this 
mean  are  relatively  large.  Alternatively,  Fig.  12(b)  presents  the  joint  probability  density  of  the  fuel  consumption 
term  away  from  the  flame  front  (at  local  temperature  T)  normalized  by  the  consumption  term  of  the  closest  point 
on  the  flame  front  {i.e.  normaly  projected  on  the  iso-surface  T  =  Tpea k).  As  can  be  observed,  the  fluctuations  are 
significantly  reduced  now.  Figure  12(c)  shows  the  same  normalized  source  term  for  the  non-unity  Lewis  number  case. 
These  plots  suggest  that  fuel  consumption  locally  scales  like  its  value  at  Tpeak,  for  both  the  unity  and  the  non-unity 
Lewis  number  flames,  i.e. 
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(b)  Le  =  1;  normalized  by  closest  point  on  (c)  Le  ^  1;  normalized  by  closest  point  on 
isoterm  T  =  Tpeak  isoterm  T  =  Tpeak 


Fig.  12.  (a) Joint  PDF  and  conditional  mean  (solid  line)  of  71-C7H16  source  term,  chc7Hi6j  normalized  by  its  peak  laminar  value  vs. 
temperature  from  the  unity  Lewis  number  simulation.  Joint  PDF  of  n- C7H16  source  term,  chc7Hi6>  normalized,  in  the  direction  normal 
to  the  isoterm  T  =  Tpeak,  by  its  value  on  this  isoterm  vs.  temperature  for  the  unity  Lewis  number  (b)  and  the  non-unity  Lewis  number 
DNS  (c).  The  non-unity  and  unity  Lewis  number  flamelet  solutions  are  also  shown  (dashed  line). 


^C7H16(T)  «  ^C7H16  (Tpeak)  T 


^C7Hi6,lam(^) 


^C7Hi6,lam(^peak) 

Therefore,  only  the  fuel  consumption  term  at  T  —  Tpeak  will  be  further  considered. 


(35) 


6.3.  Curvature  and  strain  rate 


The  mean  curvature  is  expressed  as 

n  =  V  •  n,  (36) 

where  n  =  —  VT/|VT|  is  the  normal  to  the  reaction  zone  surface,  and  is  positively  oriented  towards  the  unburnt 
gas.  The  probability  density  distribution  of  curvature  is  presented  in  Fig.  13(a)  for  both  the  unity  and  the  non-unity 
Lewis  number  cases.  First,  it  is  important  to  note  that  the  magnitude  of  curvature  is  very  large,  and  such  intense 
curvatures  are  difficult  to  observe  in  laminar  flames.  Second  and  most  importantly,  no  differential  diffusion  effect  is 
visible  in  Fig.  13(a),  suggesting  that  the  geometry  of  the  reaction  surface  is  only  influenced  by  turbulence. 


-2  -1.5  -1  -0.5  0  0.5  1  1.5  2 
Mean  curvature  x5r  H 

L7M16 

(a)  Curvature 


Strain  rate  a^^/S^ 


(b)  Strain  rate 


Fig.  13.  Probability  density  function,  on  the  reaction  surface,  of  tc,  normalized  by  the  laminar  reaction  zone  thickness  (top),  and  at, 
normalized  by  the  ratio  of  the  laminar  reaction  zone  thickness  to  the  laminar  flame  speed  (bottom). 


The  strain  rate  tangential  to  the  reaction  surface  is  computed  as,  similarly  to  Refs.  [50,51], 

at  =  V  u  —  n  •  Vw  •  n.  (37) 

where  u  is  the  velocity.  Figure  13(b)  shows  the  probability  density  of  the  normalized  strain  rate  tangential  to  the 
reaction  surface.  The  magnitude  of  the  strain  rate,  up  to  70,000  s_1,  is  very  large  and  hard  to  reproduce  in  laminar 
flames.  Again,  the  strain  rate  applied  to  both  the  unity  Lewis  number  and  the  non- unity  Lewis  number  reaction 
surface  are  similar.  Consistent  with  values  previously-reported  by  several  authors  for  methane/air  and  hydrogen/air 
flames  in  the  thin  reaction  zone  [51-54],  a  positive  mean  tangential  strain  rate  is  observed. 

With  the  information  presented  in  this  section,  it  is  interesting  to  consider  the  following  expression  (often  used  to 
predict  the  effects  of  strainrate  and  curvature) 
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^C7Hi6  + 


(38) 


£ 

^c7h16 


at5  c7h16  \ 

S°L  )' 


Under  similar  unburnt  conditions  (but  with  an  unburnt  temperature  of  353K),  C  «  1mm  [55],  which  means 
C,  / ^c7Hi6  ~  0(10).  Figure  13  shows  that  «5c7Hi6  ~  0(O.l)-0(l),and  at5c7ii16/ Sl  ~  0(1)-0(1O),  which  means 
that  £  (ft  +  at/SzJ  1  and  Eq.  38  does  not  hold  anymore. 


6.4.  Propagating  vs.  material  surface 


The  previous  results  suggest  that  differential  diffusion  virtually  does  not  affect  the  reaction  zone  geometry.  This 
is  a  consequence  of  the  high  Karlovitz  number  and  can  be  explained  by  an  analogy  to  the  work  of  Yeung  et  al.  [56] 
who  analyzed  theoretically  and  numerically  the  geometry  of  material  and  propagating  surfaces  in  the  presence  of 
homogeneous  isotropic  turbulence.  Yeung  et  al.  [56]  found  that  a  material  surface  orientates  preferentially  with  the 
strain  rate  tensor.  They  identified  a  universal  distribution  of  strain  rate  on  the  material  surface.  In  particular,  they 
observed  that  80%  of  the  strain  rate  is  positive.  They  further  argued  that  a  propagating  surface  (with  speed  w)  has 
the  same  strain  rate  distribution  as  a  material  surface,  independently  of  the  turbulent  Reynolds  number,  under  the 
condition  w/v^  <C  1,  where  vv  is  the  Kolmogorov  velocity. 

For  the  present  flames  analyzed,  81%  and  82%  of  the  strain  rate  was  found  to  be  positive  for  the  unity  and  the 
non-unity  Lewis  number  flames,  respectively.  These  results  are  consistent  with  the  ratios  S^/v^  as  listed  in  Table  3. 
These  results  are  also  consistent  with  the  series  of  DNS  of  lean  LU/air  turbulent  flames  conducted  by  Aspden  et 
al.  [5]. 


Fuel 

Ka 

Sl/vt, 

%  at  >  0 

n-CrHi6/air  (Le  =  1) 

280 

0.08 

81 

n-CrHi6/air  (Le  /  1) 

220 

0.11 

82 

H2  /  air 

28 

0.30 

80 

H2/air 

4200 

0.02 

80 

Table  3 

Sf/vv  ratios  (see  Ref.  [56])  and  curvature  and  strain  rate  statistics  for  the  flames  presented  in  this  thesis  and  for  two  lean  H2/air 
(</)  =  0.4)  flames  presented  in  [5].  vv  is  evaluated  at  Tpeak. 

The  overall  conclusion  is  that  the  reaction  zone  surface  behaves  as  a  material  surface  under  the  intense  turbulent 
field.  Even  if  the  local  consumption  rates  are  affected,  they  are  too  weak  to  alter  the  shape  of  the  flame  front.  As 
such,  the  distributions  of  curvature  and  strain  rate  are  independent  from  differential  diffusion  effects. 

6.5.  Turbulent  flame  speed 


The  turbulent  flame  speed  may  be  expressed  as 


ST 


1 

PuFc7Hi6,uF2  (h  ~  h) 


(x,y,z,t)  dVd £, 


(39) 


where  the  subscript  u  is  used  for  unburnt  conditions,  £2  ~  £1  is  the  data  collection  time.  The  values  obtained  are 
presented  in  Table  4. 

The  difference  in  the  St/ Sl  ratios  (between  unity  and  non-unity  Lewis  numbers)  highlights  an  important  differen¬ 
tial  diffusion  effect  (the  only  difference  between  the  flames).  The  unity  Lewis  number  flame  encounters  a  significantly 
larger  increase  in  flame  speed  than  the  non-unity  Lewis  number  flame,  although  both  flames  are  subjected  to  the 
same  incoming  turbulent  flow.  As  observed  by  many  authors  [57-64],  this  further  suggests  that  differential  diffusion 
has  to  be  taken  into  account  in  turbulent  flame  speed  models  (e.g.  in  the  form  of  a  Lewis  number  dependence). 

Although  the  latter  observation  is  interesting  in  itself,  a  better  understanding  of  where  this  difference  comes  from 
is  necessary  for  modeling  purposes.  Even  if  the  current  flames  are  not  in  the  (thin)  corrugated  flamelet  regime, 
the  reaction  zone  remains  thin  and  is  only  weakly  corrugated.  As  a  consequence,  a  reaction  surface  (isocontour 
T  =  Tpeak)  and  a  local  consumption  speed  can  still  be  defined.  For  the  present  analysis,  a  similar  approach  to  that 
used  by  Damkohler  for  the  corrugated  flamelet  regime  is  taken  [65]. 

As  shown  in  Fig.  12(b)  and  12(c),  the  profiles  of  cuc7Hi6  in  the  direction  normal  to  the  reaction  surface  remain 
fairly  unchanged  when  scaled  by  their  value  at  T  =  Tpeak.  Hence,  using  Eq.  35  in  Eq.  39,  the  turbulent  flame  speed 
can  be  approximated  as 


St 


^c7h16 

^C7Hi6,lam 


*  peak 


(40) 


where  (c5c7Hi6  A^c7Hi6,iam)T  k  is  surface- weighted  and  averaged  in  time.  AT  is  defined  as  the  surface  area  of  the 
reaction  zone  surface  as  shown  in  Fig.  14.  These  values  are  presented  in  Table  4  for  the  present  flames.  Once  again, 
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A 


Fig.  14.  Schematic  drawing  to  illustrate  Eq.  40.  Sp*  =  (a;c7H16/^C7Hi6,lam)  •  is  used. 

'  ’  '  ^peak 

differential  diffusion  has  a  limited  effect  on  the  turbulent  surface  area:  the  flame  acts  as  a  material  surface.  Second, 
the  values  presented  are  in  very  good  agreement  with  Eq.  40.  Third,  the  main  contribution  to  the  differences  in 
turbulent  flame  speed  due  to  differential  diffusion  appears  in  the  (cuc7Hi6/^C7Hi6,iam)T  eak  ratio. 


unity  Lei 

non-unity  Le i 

ST  (m/s) 

1.06 

0.69 

St/S°l 

3.7 

1.9 

At /A 

3.5  ±  1.0 

2.9  ±  0.8 

/  ^c7h16  \ 

\-C7H16,am/Tpeak 

1.00 

0.66 

At  /  UJc7n16  \ 

A  \  WC7H16  ,lam  /  rp 

Jpeak 

3.5 

1.9 

%/AT 

0% 

2.9% 

aioo%/ At 

56.4% 

76.6% 

Table  4 

Ratios  relevant  to  the  turbulent  flame  speed  (see  Eq.  40).  Confidence  intervals  correspond  to  plus  or  minus  one  standard  deviation. 


For  the  unity  Lewis  number  case,  only  strain  rate  has  an  effect  on  the  fuel  consumption  rate  (see  Fig.  15(a)). 
This  effect  is  approximately  symmetric  with  respect  to  the  mean  strain  rate  and  this  strain  rate  is  symmetrically 
distributed  with  respect  to  its  mean  (see  Fig.  13(b)).  Therefore,  a  close-to-unity  (c5c7Hi6/^;C7Hi6,iam)T  k  ratio  is 

expected.  This  is  consistent  with  results  previously  reported  for  turbulent  flames  with  a  close-to-unity  fuel  Lewis 
number  [50,52,66]. 

For  the  non-unity  Lewis  number  case,  the  fuel  consumption  rate  is  correlated  with  curvature  (see  Fig.  15(b)). 
While  curvature  is  symmetrically  distributed  with  respect  to  its  mean,  its  effect  on  the  fuel  consumption  term  is 
highly  non-linear  and  is  not  symmetric  with  respect  to  its  mean.  Consequently,  the  overall  effect  is  a  reduction  in 
the  average  fuel  consumption  rate,  i.e.  (uJc7u16/^c7}i16,\am)T  k  <  1. 
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Fig.  15.  Joint  probability  density  function,  on  the  isoterm  T  =  Tpeak,  of  n-C7Hi6  source  term,  c5c7h16>  normalized  by  its  peak  laminar 
value  vs.  the  normalized  tangential  strain  rate  for  Le  =  1  and  by  the  normalized  mean  curvature  for  Le  ^  1. 
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6.6.  Summary 


The  flow  chart  presented  in  Fig.  16  illustrates  the  mechanism  through  which  turbulence  and  differential  diffusion 
affect  the  overall  turbulent  flame  speed  and  summarizes  the  main  results  of  this  section. 


Turb. 
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UJp 


<^F,lam 


St 

SI 


Fig.  16.  Flow  chart  illustrating  the  mechanism  through  which  turbulence  affects  the  flame  speed  for  the  flames  presented  in  this  work. 


7.  Modeling  the  reaction  zone 


For  flames  A  and  B,  the  following  results  were  identified: 

(i)  largely  thickened  preheat  zone,  yet  thin  reaction  zone  (see  Fig.  11), 

(ii)  large  burning  rate  fluctuations  around  the  mean  profile  vs.  progress  variable, 

(iii)  for  unity  Lewis  numbers: 

<  CF\c  cjp5iam  (c) , 


(iv)  for  non-unity  Lewis  numbers: 


<Ll>f\c  > ^  UJf, lam  (c)  , 


(41) 

(42) 


(v)  for  both  cases: 


lvf  (x,  t)  7^  LdF  (c  (x,  t) )  . 


(43) 


This  section  is  directly  related  to  point  5,  i.e.  the  objective  is  to  identify  a  set  of  variables  ^  such  that  ojf  (x,  t )  can 
be  approximated  accurately  by  CoF  (c  (x,  t) ,  ^  (x,  £)),  while  keeping  the  dimensionality  of  ^  small. 


7.1.  Coordinate  transformation 

In  order  to  identify  the  set  of  variables  Vb  a  similar  approach  to  the  one  used  in  Ref.  [25]  for  non-premixed  flames 
is  taken  in  this  section.  The  following  coordinate  transformation  is  proposed: 

(x1,x2,xs,t)  ->►  (c(x1,x2,xs,t)  ,c2  {x\ , x2 , x3 ,  t) ,  c3  (xx , x2 , x3 ,  t) , t)  ,  (44) 

with 

Vc  •  Vc2  =  0,  and  Vc  •  VC3  =  0,  (45) 

i.e.  the  variables  C2  and  C3  lie  in  the  surface  of  constant  c.  The  only  assumption  made  in  this  section  is  that  the 
Jacobian  of  the  transformation  is  not  singular.  Applying  the  coordinate  transformation  to  Eq.  5,  the  following 
transformed  equations  for  the  species  mass  fractions  are  obtained: 
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where  y  =  2<a|Vc|2  and  =  2ct|Vc/e|2  for  k  =  2, 3  are  the  dissipation  rates  of  c,  C2,  and  C3,  respectively.  In  the  above 
equations,  the  Lagrangian  transport  corresponds  to  the  convection  of  the  scalar,  in  the  C2-C3-direetion,  i.e.  in  the 
surface  of  iso-c,  induced  by  the  Lagrangian  transport  of  these  C2  and  C3  coordinates.  The  normal  terms  correspond 
to  convection  or  diffusion  of  the  scalar  in  the  c-direction,  whereas  the  tangential  terms  correspond  to  convection  or 
diffusion  of  the  scalar  in  the  C2-C3-direetion. 


7.2.  Flamelet  equations  -  Unity  Lewis  number 


Recall  that,  in  the  regime  considered,  the  reaction  zone  is  considered  to  be  thin.  It  is  therefore  assumed  that,  in 
the  vicinity  of  the  reaction  zone,  d/dck  «  d/dc,  for  k  =  2,3.  Furthermore,  assuming  the  dependence  in  r  in  the 
transformed  coordinate  system  is  negligible,  the  flamelet  equations  are  obtained: 


convection 


PX  d2!) 
'  2  dc2 

diffusion 


chemical  source 


(47) 


These  equations  form  a  system  of  ordinary  differential  equations  in  which  the  Y;L  are  the  unknowns,  c  is  the  variable, 
and  the  dissipation  rate  y(c)  is  the  parameter.  Each  of  the  terms  in  the  above  equations  are  only  a  function  of  c  and 
X- 

Given  the  form  of  the  flamelet  equations,  the  fuel  burning  rate  should  solely  be  a  function  of  the  progress  variable 
c  and  its  dissipation  rate  y.  In  contrast  with  Fig.  15(a)  and  15(b)  shown  previously,  Fig.  17  presents  the  joint 
probability  density  function  of  the  fuel  consumption  rate  vs.  the  dissipation  rate  of  the  progress  variable.  It  is  clear 
that  the  fuel  burning  rate  is  far  more  correlated  with  x  than  with  strain-rate.  Physically,  this  means  that  turbulence 
affects  the  fuel  burning  rate  by  compressing  or  extending  the  isosurfaces  of  the  progress  variable.  The  figure  also 
shows  that  the  solution  of  the  flamelet  equations  (Eq.  47)  provide  a  very  good  estimate  of  the  dependence  of  the 
chemical  source  term  on  the  dissipation  rate. 


7.3.  Flamelet  equations  -  Non-unity  Lewis  number 

Once  again,  we  make  the  same  assumptions  of  local  one-dimensionality  in  phase  space,  i.e.  d/dck  «  d/dc ,  for 
k  =  2,3,  and  steady  state  in  phase  space,  i.e.  d/dr  =  0,  but  no  assumption  is  made  on  the  species  Lewis  numbers. 
Let  £  =  V  •  (paVc).  Then,  non-unity  Lewis  number  flamelet  equations  become 
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2.5 


Fig.  17.  Comparison  between  the  fuel  burning  rate  predicted  by  Cjp  fgm  (c,  x)  and  <  |  (c,  x)  >  as  a  function  of  the  dissipation  rate, 

normalized  by  their  peak  laminar  values,  given  c  ==  cpeak-  The  joint  PDF  on  the  isosurface  of  c  =  cpeak  for  the  unity  Lewis  number  flame 
B  is  also  shown  for  comparison. 
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(48) 


In  this  set  of  equations,  a  second  parameter,  £  (c),  is  involved.  It  corresponds  to  the  diffusion  of  c  in  the  original 
coordinate  system.  Note  that,  as  one  would  expect,  this  parameter  includes  curvature  effects. 

Figure  18  shows  the  probability  density  function  of  the  point-wise  comparison  in  flame  B  between  the  predicted 
and  the  actual  fuel  burning  rate,  for  three  different  optimal  estimators:  <  cup|c  >,  <  cup|c,  x  >?  and  <  cup|c,  %,£  >. 
The  inclusion  of  the  third  parameter  leads  to  obvious  improvements. 


Jam, peak 

(a)  <  Cjf\c  >  vs.  t) 


lam,  peak 

(b)  <  ujf \c,x  >  vs.  <hp(x,  t) 


Jam,  peak 

(c)  <  d;F|c,x,£  >  VS.  vF(x,t) 


Fig.  18.  Joint  PDF  of  the  comparison  between  the  predicted  fuel  burning  rate  and  the  actual  burning  rate  in  the  higher- Ka  non-unity 
Lewis  number  flame. 


7.4.  Summary 

While  the  flames  considered  in  this  work  have  a  thin  reaction  zone,  large  fluctuations  around  the  mean  burning  rate 
(conditional  on  c)  were  observed.  In  order  to  model  these  fluctuations,  a  coordinate  transformation  was  performed, 
and  the  following  results  were  identified: 

(i)  for  the  unity  Lewis  number  flames, 

(a)  the  turbulent  flames  can  be  well  represented  by  a  set  of  one-dimensional  (in  c- space)  flamelet  equations 
parameterized  by  x  (c)? 
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(b)  the  fuel  burning  rate  uop  (x,  t)  can  be  approximated  accurately  by  ujf  (c  (x,  t) ,  x  (x,  £)), 

(c)  the  fuel  burning  rate  is  well  predicted  a  priori  by  the  solution  to  the  set  of  flamelet  equations, 

(ii)  for  the  non-unity  Lewis  number  flames, 

(a)  the  fuel  burning  rate  Cop  (x,  t)  can  be  approximated  accurately  (only  to  some  extent  in  the  higher-Pa 
flame)  by  uF  (c  (x,  t) ,  x  (x,  t) ,  £  (x,  t)), 

(b)  the  fuel  burning  rate  cannot  be  predicted  a  priori  by  the  solution  to  the  set  of  flamelet  equations,  unless 
a  model  for  the  effective  diffusivities  as  a  function  of  the  progress  variable  is  considered  (future  work), 

(c)  the  source  term  of  the  progress  variable  (relevant  to  a  turbulent  flame  simulation  with  this  reduced 
chemistry)  can  be  closed  by  an  a  posteriori  DNS-generated  manifold. 

Finally,  it  is  important  to  note  that  the  assumption  of  a  thin  reaction  zone  was  made  to  obtain  the  premixed 
flamelet  equations.  While  this  assumption  is  valid  in  the  thin  reaction  zone  regime,  it  may  not  be  the  case  for  higher 
turbulence  intensities  corresponding  to  the  distributed  burning  regime. 


8.  Overview  of  turbulent  flow 


8.1.  Conditional  averaging 


Turbulence  quantities  are  expected  to  vary  through  the  flame  based  upon  the  local  thermodynamic  properties  of 
the  fluid  (such  as  density  and  viscosity).  In  a  curved  and  instantaneously  transient  flame,  these  quantities  correlate 
better  with  a  flame  progress  variable,  denoted  as  (7,  compared  to  the  spatial  coordinate  x.  For  this  reason,  averages 
are  conditioned  on  C  and  denoted  as 

W  C)  (49) 

The  form  of  the  progress  variable  chosen  is  C  =  Yh2o  +  Yh2  +  Yco2  +  as  it  tracks  the  flame  evolution  through 

the  preheat  and  reaction  zones.  The  progress  variable  range  is  standardized  by  considering  C  =  C/Cmaa,  so  that 
0  represents  the  reactants  and  1  represents  the  products.  By  use  of  this  conditional  averaging,  we  define  the  local 
Kolmogorov  time,  length,  and  velocity  scales  as 


Tr,{C)  = 
V(C)  = 

UV(C)  = 


m)  ' 

(He)5 


1/4 


((e\C)(v\C))1/A, 


(50) 

(51) 

(52) 


the  local  Karlovitz  number  as  Ka(C )  =  Tp/r^C),  and  the  dissipation  rate  conditioned  on  C  as  (e\C)  =  (r'  : 
S"\C) / (p\C).  Additionally,  we  define  a  quantity  involving  the  flame  density  change,  7 (C)  =  Ap/(p|C),  where  A p  = 
Pu  —  Pb,  which  will  be  used  in  the  subsequent  analysis. 


8.2.  Vorticity  overview 

First,  the  transformation  of  vorticity  is  qualitatively  observed  by  plotting  the  vorticity  magnitude  through  the 
flame.  Figure  19  displays  instantaneous  2D  contours  of  cases  A,  B,  C,  D  and  E.  The  change  in  vorticity  is  dramatic 
as  the  magnitude  is  greatly  suppressed  through  the  flame.  Preliminary  observation  suggests  the  vorticity  is  reduced 
to  a  lesser  extent  in  cases  C  and  D,  which  have  a  similar  Kau  as  B  and  E,  respectively,  but  a  higher  unburnt 
temperature. 


9.  Enstrophy  budget 

The  enstrophy,  uj2  =  u?  •  Cc?,  transport  equation  is  derived  from  the  momentum  equation  as 

-V  •  T  ]  +  Cc?  •  V  X  — , 
p  )  p 

where  D / Dt  is  the  material  or  total  derivative.  Each  term  on  the  right  hand  side  is  associated  with  a  specific 
physical  processes:  vortex  stretching,  dilatation,  baroclinic  torque,  viscous  dissipation,  and  forcing,  respectively. 
Vortex  stretching,  viscous  dissipation,  and  forcing  are  active  in  constant  density  flows,  while  dilatation  and  baroclinic 
torque  arise  here  only  due  to  the  presence  of  the  flame.  The  change  of  fluid  properties  (such  as  density  and  viscosity) 
within  a  premixed  flame  alters  the  enstrophy  of  the  incoming  turbulence  through  these  five  terms. 

In  the  following,  a  scaling  for  each  term  is  provided.  Only  the  first  term  is  detailed. 


1  Duj2 


2  Dt 


=  UJ  •  S  •  Cc? 


LO 


•  uo2  (V  •  u)  +  •  (Vp  x  VP)  +  u?  •  V 
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Fig.  19.  Two-dimensional  slices  of  vorticity  magnitude  for  cases  A,  B,  C,  D  and  E.  Each  figure  corresponds  to  a  region  of  size  L  x  5 L 
and  the  flow  direction  is  from  bottom  to  top.  The  contours  range  between  [0;  2.1e5  s_1],  [0;  6.2e5  s_1],  [0;  5.5e6  s_1],  [0;  2.1e7s-1],  and 
[0;  1.7e6  s-1l,  respectively.  Blue  and  red  contours  represent  the  extent  of  the  turbulent  flame  brush  defined  as  the  iso-surfaces  of  C  =  0.05 
and  C  =  0.8. 


9.1.  Vortex  stretching 

In  dimensional  form,  the  vortex  stretching  term  can  vary  by  more  than  an  order  of  magnitude  across  the  flame  and 
varies  by  several  orders  of  magnitude  between  the  different  cases  (Fig.  20a).  Scaling  of  this  term  requires  estimates 
for  vorticity  and  the  rate  of  strain  tensor. 

In  the  present  configuration  of  high  Karlovitz  number  ( Kau  >70),  we  estimate  the  magnitude  of  S  with  1/r^. 
This  result  is  consistent  with  homogeneous,  isotropic  turbulence  where  S'  :  S'  also  scales  as  (l/r2)  [67].  As  a  spatial 
gradient  of  velocity,  like  the  rate  of  strain  tensor,  vorticity  is  estimated  as  1/r^  and  enstrophy  with  1/r2.  Previous 
experimental  and  numerical  work  supports  a  correlation  of  vorticity  with  the  Kolmogorov  time  scale  under  conditions 
of  homogeneous,  isotropic  turbulence  [68-70].  The  above  analysis  results  in  the  scaling, 

u>  •  S  •  u)  oc  (53) 

Normalized  according  to  the  above  expression,  the  vortex  stretching  terms  for  each  case  collapse  to  a  fairly 
constant  value  close  to  0.15,  which  is  the  same  value  obtained  in  DNS  of  homogeneous,  isotropic  turbulence  (Fig. 
20b).  The  success  of  the  normalization  suggests  that  within  the  flame  the  vortex  stretching  term  behaves  similar  to 
homogeneous,  isotropic  turbulence  in  the  limit  of  high  Karlovitz  number. 


9.2.  Dilatation 


Dilatation  in  the  present  configuration  is  due  only  to  the  effects  of  the  flame.  Scaling  the  dilatation  term  requires 
estimates  of  enstrophy  and  the  divergence  of  velocity.  This  leads  to  the  following  scaling, 


cj2(V  •  u)  oc  7 


(54) 
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(a)  Dimensional. 


(b)  Normalized. 


Fig.  20.  Vortex  stretching  term  in  dimensional  and  normalized  form.  Oscillations  at  the  smallest /largest  values  of  C  are  due  to  the  region 
where  turbulence  forcing  is  inactive.  The  line  Re\  is  calculated  from  the  previous  simulations  of  homogeneous,  isotropic  turbulence 
by  [44]. 


9.3.  Baroclinic  torque 

Baroclinic  torque,  like  dilatation,  is  only  present  due  to  the  density  variations  within  the  flame.  Likewise,  it  tends 
toward  zero  in  the  reactants  and  products.  The  following  scaling  for  baroclinic  torque  is  proposed, 

~2  '  (Vp  x  VP)  oc  (55) 

P 2  l  ft 

When  normalized  according  to  the  above  expression,  the  variation  in  the  peak  value  between  the  six  cases  reduces 
from  five  orders  of  magnitude  to  a  factor  of  2. 


9.4.  Viscous  dissipation 

Viscous  dissipation,  like  vortex  stretching,  is  present  in  constant  density  flow.  As  with  several  quantities  discussed 
prior,  in  the  case  of  high  iLa,  viscous  dissipation  scales  with  the  turbulent  quantities,  specifically, 

“■vxGv4V  <56) 

In  the  the  high  Reynolds  number  limit  of  homogeneous,  isotropic  turbulence,  [71]  propose  a  similar  scaling. 


9.5.  Normalized  enstrophy  transport  equation 

The  above  scaling  estimates  are  used  to  propose  a  normalization  of  the  entire  enstrophy  transport  equation,  which 
is  then  given  by, 


2  ~  KaTi  4-  7^2  +  a^fV KaT^  +  KaT \  +  ^Da  (^7) 

where,  t  =  t/rp,  a)2  =  (cu2|C)r2,  and  the  Damkohler  number  is  Da  =  t0/tf.  As  written  above,  T\  is  vortex 

stretching,  T2  is  dilatation,  T3  is  baroclinic  torque,  T4  is  viscous  dissipation,  and  T5  is  the  forcing  with  each  component 
normalized  according  to  the  above  discussion.  The  preceding  analysis  in  this  section  shows  that  the  normalized  terms 
obtain  nearly  constant  values  (or  a  constant  peak  value)  through  the  flame  and  across  conditions  with  T\  ~  0.15, 

0  <  T2  <  0.7,  0  <  T3  <  0.1,  —0.25  <  T4  <  —0.15,  and  T5  ~  1,  at  high  Kau.  The  scaling  in  Eq.  57  suggests  that 
as  the  Karlovitz  number  increases,  vortex  stretching  and  dissipation  increase  in  magnitude  relative  to  baroclinic 
torque,  dilatation,  and  forcing.  This  explains  the  observed  relative  decrease  of  baroclinic  torque  and  dilatation  as 
Ka  increased  in  the  results  of  Hamlington  et  al.  [72]. 
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c 

(a)  Dimensional. 


(b)  Normalized. 


Fig.  21.  Enstrophy  in  dimensional  and  normalized  form.  In  normalized  form,  enstrophy  has  less  variation  through  the  flame  as  the 
Karlovitz  number  increases.  The  line  Re^  is  calculated  from  the  previous  simulations  of  homogeneous,  isotropic  turbulence  by [44]. 


9.6.  Kolomogorov’s  first  hypothesis 


In  the  limit  of  high  Ka ,  enstrophy  transport  results  in  a  local  balance  between  production  and  dissipation. 
This  is  also  the  case  for  homogeneous,  isotropic  turbulence[73,71].  This  implies  that  enstrophy  should  scale  with 
the  Kolmogorov  time  scale,  which  is  confirmed  in  Fig.  21,  with  a  proportionality  constant  close  to  unity.  From 
homogeneous,  isotropic  turbulence  it  can  be  shown  that 


L d2  =  15 


2 


(58) 


so  that  a  normalized  value  of  unity  is  expected  [67].  This  demonstrates  that  it  is  the  change  in  kinematic  viscosity, 
as  opposed  to  the  effect  of  the  density  change  through  dilatation  and  baroclinic  torque,  which  drives  the  enstrophy 
transformation  through  the  flame. 

These  results  shed  light  on  the  validity  of  Kolmogorov’s  first  similarity  hypothesis  within  premixed  flames,  since 
enstrophy  is  characteristic  of  the  smallest  turbulent  scales.  Though  the  fluid  properties  and  turbulence  characteristics 
vary  widely  across  the  flame  and  between  the  present  cases,  enstrophy  is  found  to  vary  only  as  a  function  of  r„,  or 
equivalently  v  and  e  alone.  This  supports  that  Kolmogorov’s  first  similarity  hypothesis  is  valid  for  sufficiently  high 
Karlovitz  number  premixed  flames. 


10.  Impact  of  models  and  assumptions  on  turbulence  evolution 

The  above  conclusions  are  obtained  using  finite-rate  chemistry  and  non- unity  Lewis  numbers.  Turbulent  combus¬ 
tion  research  often  employs  other  chemical  and  transport  models. 

10.1.  Effects  of  unity- Lewis  number  assumption 

To  test  the  effects  of  transport  models,  case  B  is  repeated  setting  all  Lewis  numbers  to  unity,  referred  to  as  Bi. 
Figure  22a  and  b  show  that  varying  the  transport  model  has  negligible  effects  on  the  dimensional  and  normalized 
mean  enstrophy.  As  vortex  stretching  and  viscous  dissipation  dominate  the  transport  of  enstrophy  at  high  Karlovitz 
number,  the  mean  enstrophy  evolves  similarly  in  the  two  cases.  It  is  important  to  note  that  this  analysis  is  in  the 
absence  of  thermo-diffusive  instabilities,  and  these  results  may  not  extend  to  flames,  such  as  lean  hydrogen/air, 
where  these  instabilities  occur. 

10.2.  Effect  of  chemical  models 

To  investigate  the  implications  of  chemical  models,  case  Bi  is  repeated  using  two  alternative  chemical  mechanisms: 
one-step  (Bos,i)  and  tabulated  chemistry  (Bxab,i)-  In  order  to  focus  on  the  effects  of  the  chemical  model,  unity  Lewis 
numbers  are  used  in  all  three  simulations  in  order  to  eliminate  effects  due  differences  in  the  transport  models. 

The  results  obtained  with  each  model  are  compared  primarily  through  the  mean  enstrophy.  Figure  23a  and  b 
shows  that  the  chemical  models  induce  little  differences  in  the  dimensional  and  normalized  mean  enstrophy.  The 
agreement  observed  in  the  mean  enstrophy  further  emphasizes  that  the  smallest  turbulent  scales  evolve  through  the 
flame  with  the  local  fluid  properties,  independent  of  the  tested  transport  or  chemical  models. 
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(a)  Dimensional  enstrophy.  (b)  Normalized  enstrophy. 

Fig.  22.  Comparison  of  case  B  with  constant  non-unity (B)  and  unity  Lewis  numbers(Bi). 


(a)  Dimensional  enstrophy.  (b)  Normalized  enstrophy. 


Fig.  23.  Comparison  of  case  Bi  with  finite-rate,  tabulated,  and  one-step  chemistry. 


10.3.  Extension  to  larger  Reynolds  numbers 


While  this  study  considers  a  wide  range  of  Karlovitz  numbers,  the  integral  length  scale  in  each  DNS  discussed  thus 
far  is  limited  to  approximately  the  flame  thickness.  It  is  important  to  consider  whether  these  results  are  applicable 
to  larger  values  of  Iq/If ,  or  equivalently,  larger  Reynolds  numbers. 

The  portion  of  enstrophy  contained  in  scales  smaller  than  the  flame  thickness  may  be  determined  by  considering 
the  model  spectrum  of  Pope  for  the  high  Re  limit  [67], 


E(k)  =  1.5  e2/3/t“5/3 


K/ln 


[(nl0)2  +  6.78]1/2 


11/3 


exp 


-5.2 


([M4  +  0.44]1/4-0.4)J, 


(59) 


and  evaluating  the  vorticity  spectrum  as  Q(k)  =  2 k2E(k).  Considering  an  infinite  Re,  at  least  80%  of  the  enstrophy  is 
contained  in  scales  smaller  than  the  flame  if  If/ V  >  40.  This  is  approximately  the  condition  of  case  B.  Said  otherwise, 
if  case  B  had  an  arbitrarily  large  Re  but  the  same  Kau,  then  only  20%  of  the  enstrophy  would  be  in  scales  larger 
than  the  flame  thickness.  In  summary,  any  high  Reynolds  number  turbulent  flame  with  the  same  Karlovitz  number 
of  the  present  cases  would  have  nearly  the  same  fraction  of  vorticity  contained  in  scales  smaller  than  the  flame.  The 
present  results  should  thus  remain  valid  for  larger  integral  to  flame  length  scale  ratios. 

An  additional  DNS,  labeled  case  B^abl,  is  performed  to  verify  if  the  results  are  indeed  independent  of  the 
Reynolds  number.  This  simulation  has  the  exact  same  conditions  as  case  BTab,b  but  a  higher  Reynolds  number, 
Ret  =  1150.  The  increase  in  Ret  is  accomplished  by  increasing  both  L  as  well  as  id,  in  order  to  maintain  the  same 
Kau  number,  resulting  in  a  four-fold  increase  in  the  integral  length  scale  ( Iq/If  =  4  and  u' /Sl  =  33).  Figure  24 
presents  instantaneous  density  contour  plots  from  BTab,i  and  B^ab  v  Figure  25  shows  that  increasing  the  turbulent 
Reynolds  number  introduces  only  small  changes  to  these  quantities. 
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(b)  case  BTabl 


Fig.  24.  Two-dimensional  slices  of  density  from  cases  Bxab,i  and  B^ab  Each  figure  represents  a  region  of  size  L  x  6 L  and  both  are 
scaled  to  match  physical  length  scales. 


(a)  Dimensional  enstrophy. 


(b)  Normalized  enstrophy. 


Fig.  25.  Comparison  of  case  B^ab,!  with  the  case  of  a  larger  Reynolds  number,  B^ab  1 
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